#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat May  8 03:08:38 2021

@author: sarupa
"""
#Imports 
import pandas as pd
import copy
import numpy as np
import tensorflow as tf
from tensorflow import keras
from matplotlib import pyplot as plt
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense, Dropout, SimpleRNN,Flatten,TimeDistributed
from sklearn.preprocessing import MinMaxScaler
from LSTM_model import ObtainlstmModel
from sklearn.metrics import mean_squared_error
from casadi import *
from numpy import linalg as la
from control import *
import timeit
import time
from numpy import random
from random import seed

start = timeit.default_timer()

# Instruction: while running with saved model comment the lines from 122 to 227
# And comment the lines from 236 to 249

# Define RMSE
def rmse(predictions, targets):

    differences = (predictions - targets)/targets                      

    differences_squared = differences ** 2                    

    mean_of_differences_squared = differences_squared.mean()  

    rmse_val = np.sqrt(mean_of_differences_squared)*100         

    return rmse_val 

# Load the data
df=pd.read_csv("Final_data_set.csv") # This data has already been resampled
# Resample the data
dataset_in=df.iloc[:,].values
dataset=np.zeros((len(dataset_in)-1,11))
dataset[:,0]=dataset_in[1:,7]
dataset[:,1]=dataset_in[1:,2]
dataset[:,2]=dataset_in[1:,5]
dataset[:,3]=dataset_in[1:,8]
dataset[:,4:11]=dataset_in[:-1,9:16]

# The min-max scaling approach
sc= MinMaxScaler(feature_range=(0,1))
dataset_scaled = sc.fit_transform(dataset)

#Extract the scaling parameters
data_min = sc.data_min_
data_max = sc.data_max_

# Scale the test data using the scaling parameters of the original data
# dataset_new_scaled = (dataset_new_reduced - data_min)/(data_max - data_min)

# Split the original data into training, validation, and testing data sets
data_train = dataset_scaled[:200000,:]
data_validation  = dataset_scaled[400000:475000,:]
data_test = dataset_scaled[475000:480000,:]

# Prepare the data for the RNN
sequenceLength = 2
numberOfFeatures = 11
numberOfOutputs = 1
sequenceLength_out=2
number_of_states= 1

# Prepare the training data
X_train, y_train = [], []

for i in range(sequenceLength, len(data_train)-sequenceLength):
    X_train.append(data_train[i-sequenceLength:i,0:numberOfFeatures])
    y_train.append(data_train[i:i+numberOfOutputs,0:number_of_states])

X_train, y_train = np.array(X_train), np.array(y_train)

X_train = X_train.reshape(X_train.shape[0], X_train.shape[1], numberOfFeatures)
y_train = y_train.reshape(y_train.shape[0], y_train.shape[1], numberOfOutputs*number_of_states)

# Prepare the validation data
X_val, y_val = [], []

for i in range(sequenceLength, len(data_validation)-sequenceLength):
    
    X_val.append(data_validation[i-sequenceLength:i,0:numberOfFeatures])
    y_val.append(data_validation[i:i+numberOfOutputs,0:number_of_states])

X_val, y_val = np.array(X_val), np.array(y_val) 

X_val = X_val.reshape(X_val.shape[0], X_val.shape[1], numberOfFeatures)
y_val = y_val.reshape(y_val.shape[0], y_val.shape[1], numberOfOutputs*number_of_states)


# Prepare the test data for testing of the LSTM model

X_test, y_test = [], []

for i in range(sequenceLength, len(data_test)-sequenceLength):
    
    X_test.append(data_test[i-sequenceLength:i,0:numberOfFeatures])
    y_test.append(data_test[i:i+numberOfOutputs,0:number_of_states])

X_test, y_test = np.array(X_test), np.array(y_test)

X_test = X_test.reshape(X_test.shape[0], X_test.shape[1], numberOfFeatures)
y_test = y_test.reshape(y_test.shape[0], y_test.shape[1], numberOfOutputs*number_of_states)

# Building the LSTM model 

model=Sequential()
model.add(LSTM(units=50, input_shape=(sequenceLength,numberOfFeatures), return_sequences=True))
model.add(LSTM(units=50, return_sequences=True))
model.add(TimeDistributed(Dense(numberOfOutputs*number_of_states)))
model.compile(loss='mse', optimizer='adam')#metrics='mean_absolute_error')

history = model.fit(X_train, y_train, epochs=10 , batch_size=100, validation_data=(X_val, y_val))#, callbacks=[callback])#callbacks=[callbacks])

# Prediction of test data
predicted_trajectory=model.predict(X_test)
predicted_trajectory=predicted_trajectory[:,-1,:]
y_test = y_test.reshape(y_test.shape[0],-1)
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                
# Rescale to original data
y_test_unscaled = y_test*(data_max[0:number_of_states]-data_min[0:number_of_states]) + data_min[0:number_of_states]
predicted_trajectory_unscaled = predicted_trajectory*(data_max[0:number_of_states]-data_min[0:number_of_states]) + data_min[0:number_of_states]

# Data plot 
plt.figure(1)
str_output = ['CB3']
p1=1
for i in range(0): 
    
    plt.subplot(p1,1,i+1)
    plt.plot(predicted_trajectory_unscaled[0:len(data_test),0],'b')
    plt.plot(y_test_unscaled[0:len(data_test),0],'m')
    plt.ylabel("Tempstep " + str(i+1))
    plt.ylabel(str_output[i])
    plt.legend(['Data','RNN'])
    plt.grid()
    plt.xlabel("Timestep")
    if i == 0:
        plt.title('PredictedTrajectory_soft sensor Vs Actual trajectory')
    
rmse_ekf=rmse(predicted_trajectory_unscaled, y_test_unscaled)
print('rmse_softsensor',rmse_ekf)
print("MSE of model:%.5f"% mean_squared_error(predicted_trajectory_unscaled, y_test_unscaled))


# Extract the weights and obtain the mathematical expression of the LSTM [ START ]---------------------------------------------------

#Layer 1
#The number of units
units_1 = int(model.layers[0].trainable_weights[0].shape[1]/4)

W_1=model.layers[0].get_weights()[0] # The Kernel, the weighting matrix for the input
U_1=model.layers[0].get_weights()[1] # The recurrent kernel, the weighting matrix for the hidden states
b_1=model.layers[0].get_weights()[2] # The bias matrix

#Obtain the matrices for the gates and the cell state

#The kernel matrices for the gates and the cell state
W_i_1 = W_1[:,:units_1] # The input gate
W_f_1 = W_1[:,units_1: units_1*2] # The forget gate
W_c_1 = W_1[:, units_1*2:units_1*3] # The cell state
W_o_1 = W_1[:, units_1*3:] # The output gate

#The recurrent kernel matrices for the gates and the cell state

U_i_1 = U_1[:,:units_1] # The input gate
U_f_1 = U_1[:,units_1:units_1*2] # The forget gate
U_c_1 = U_1[:, units_1*2:units_1*3] # The cell state
U_o_1 = U_1[:, units_1*3:] # The output gate


# The bias matrices for the gates and the cell state

b_i_1 = b_1[:units_1] # The input gate
b_f_1 = b_1[units_1:units_1*2] # The forget gate
b_c_1 = b_1[units_1*2:units_1*3] #The cell state
b_o_1 = b_1[units_1*3:] # The output gate

#Layer 2
#The number of units
units_2 = int(model.layers[1].trainable_weights[0].shape[1]/4)

W_2=model.layers[1].get_weights()[0] # The Kernel, the weighting matrix for the input
U_2=model.layers[1].get_weights()[1] # The recurrent kernel, the weighting matrix for the hidden states
b_2=model.layers[1].get_weights()[2] # The bias matrix

#Obtain the matrices for the gates and the cell state

#The kernel matrices for the gates and the cell state
W_i_2 = W_2[:,:units_2] # The input gate
W_f_2 = W_2[:,units_2: units_2*2] # The forget gate
W_c_2 = W_2[:, units_2*2:units_2*3] # The cell state
W_o_2 = W_2[:, units_2*3:] # The output gate

#The recurrent kernel matrices for the gates and the cell state

U_i_2 = U_2[:,:units_2] # The input gate
U_f_2 = U_2[:,units_2:units_2*2] # The forget gate
U_c_2 = U_2[:, units_2*2:units_2*3] # The cell state
U_o_2 = U_2[:, units_2*3:] # The output gate


# The bias matrices for the gates and the cell state

b_i_2 = b_2[:units_2] # The input gate
b_f_2 = b_2[units_2:units_2*2] # The forget gate
b_c_2 = b_2[units_2*2:units_2*3] #The cell state
b_o_2 = b_2[units_2*3:] # The output gate

# The Dense layer
V_m=model.layers[2].get_weights()[0]
b_m=model.layers[2].get_weights()[1]


# Extract the weights and obtain the mathematical expression of the LSTM [ end ]---------------------------------------------------


# Import saved model: while running with saved model comment the lines from 122 to 227
# And comment the lines from 236 to 249
# Load the weights, units and LSTM properties
# W_1 = np.loadtxt('./Saved Models ss/W_1.txt')
# W_2 = np.loadtxt('./Saved Models ss/W_2.txt')
# U_1 = np.loadtxt('./Saved Models ss/U_1.txt')
# U_2 = np.loadtxt('./Saved Models ss/U_2.txt')
# b_1 = np.loadtxt('./Saved Models ss/b_1.txt')
# b_2 = np.loadtxt('./Saved Models ss/b_2.txt')
# V_m = np.loadtxt('./Saved Models ss/V_m.txt')
# b_m = np.loadtxt('./Saved Models ss/b_m.txt')
# data_max=np.loadtxt('./Saved Models ss/data_max.txt')
# data_min=np.loadtxt('./Saved Models ss/data_min.txt')
# y_test_unscaled = np.loadtxt("./Saved Models ss/y_test_unscaled.txt")
# y_test_unscaled= np.reshape(y_test_unscaled, (len(y_test_unscaled),1))
# units_1=50
# units_2=50

# #performance check

np.random.seed(5) 
sigma_v_m=0.005
Nsim = y_test.shape[0]
v = sigma_v_m*random.randn(Nsim,sequenceLength,3)
X_test1=copy.deepcopy(X_test)
X_test1[:,:,1:4]= X_test[:,:,1:4]+v
X_test1[0,:,0]= X_test[0,:,0]*0.8

# The mathematical expression for the LSTM and then compare the results with the model.predict 
y_rnn_model= np.zeros((len(y_test_unscaled),1))
x = X_test1[0,:,0:number_of_states]
u = X_test1[0,:,number_of_states:numberOfFeatures]
y_rnn_model[0:sequenceLength]=x
start2 = time.perf_counter()
for k in range(sequenceLength,y_test.shape[0]-sequenceLength):
    print(k)
    y_rnn_model[k,:] = ObtainlstmModel(W_1, U_1, b_1, W_2, U_2, b_2, V_m, b_m, units_1, units_2, sequenceLength, x,u)
    x=y_rnn_model[k-1:k+1,:] 
    u = X_test1[k-1,:,number_of_states:numberOfFeatures]
y_rnn_model_unscaled = y_rnn_model*(data_max[0:number_of_states]-data_min[0:number_of_states]) + data_min[0:number_of_states]
finish2 = time.perf_counter()
print('\n Finished online in time',finish2 - start2,'seconds')
plt.figure(2)
str_output = ['CB3']
p1=1
for i in range(p1): 
    
    plt.subplot(p1,1,i+1)
    plt.plot(y_test_unscaled[:400],'b')
    plt.plot(y_rnn_model_unscaled[:400],'r')
    plt.ylabel("Tempstep " + str(i+1))
    plt.ylabel(str_output[i])
    plt.legend(['DATA', 'RNN'])
    plt.grid()
    plt.xlabel("Timestep")
    if i == 0:
        plt.title('PredictedTrajectory_Rnn Vs Actual trajectory')


# To save the model parameters
# np.savetxt('./Saved Models ss/W_1.txt', W_1)
# np.savetxt('./Saved Models ss/W_2.txt', W_2)
# np.savetxt('./Saved Models ss/U_1.txt', U_1)
# np.savetxt('./Saved Models ss/U_2.txt', U_2)
# np.savetxt('./Saved Models ss/b_1.txt', b_1)
# np.savetxt('./Saved Models ss/b_2.txt', b_2)
# np.savetxt('./Saved Models ss/V_m.txt', V_m)
# np.savetxt('./Saved Models ss/b_m.txt', b_m)
# np.savetxt('./Saved Models ss/data_min.txt', data_min)
# np.savetxt('./Saved Models ss/data_max.txt', data_max)
# np.savetxt('./Saved Models ss/predicted_trajectory_unscaled.txt', predicted_trajectory_unscaled)
# np.savetxt('./Saved Models ss/y_rnn_model_unscaled.txt', y_rnn_model_unscaled)
# np.savetxt('./Saved Models ss/y_test_unscaled.txt', y_test_unscaled)
# np.savetxt("./Saved Models ss/y_test.txt", y_test)
# np.savetxt("./Saved Models ss/y_test_unscaled.txt", y_test_unscaled)
# np.savetxt("./Saved Models ss/X_test.txt", X_test.reshape(X_test.shape[0], int(X_test.shape[1]*X_test.shape[2])))


# The RMSE is reported for different initial guess and different random seed taking average
# This is a sample calculation

rmse_ekf=rmse(y_rnn_model_unscaled[0:len(y_test_unscaled)-5,:],y_test_unscaled[0:len(y_test_unscaled)-5])
print('rmse_ekf',rmse_ekf)

finish = time.perf_counter()
print('\n Finished in time',finish - start,'seconds')


  
